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A full energy and flavor-dependent analysis of the three-year high-energy IceCube neutrino events 
is presented. By means of multidimensional fits, we derive the current preferred values of the high- 
energy neutrino flavor ratios, the normalization and spectral index of the astrophysical fluxes, and 
the expected atmospheric background events, including a prompt component. A crucial assumption 
resides on the choice of the energy interval used for the analyses, which significantly biases the 
results. When restricting ourselves to the ~30 TeV—3 PeV energy range, which contains all the 
observed IceCube events, we find that the inclusion of the spectral information improves the fit 
to the canonical flavor composition at Earth, (1:1: 1)®, with respect to a single-energy bin 
analysis. Increasing both the minimum and the maximum deposited energies has dramatic effects 
on the reconstructed flavor ratios as well as on the spectral index. Imposing a higher threshold of 
60 TeV yields a slightly harder spectrum by allowing a larger muon neutrino component, since above 
this energy most atmospheric tracklike events are effectively removed. Extending the high-energy 
cutoff to fully cover the Glashow resonance region leads to a softer spectrum and a preference for 
tan neutrino dominance, as none of the expected electron antineutrino induced showers have been 
observed so far. The lack of showers at energies above 2 PeV may point to a broken power-law 
neutrino spectrum. Future data may confirm or falsify whether or not the recently discovered 
high-energy neutrino fluxes and the long-standing detected cosmic rays have a common origin. 
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I. INTRODUCTION 

The observation of high-energy neutrinos in the IceCube detector at the South Pole [1-3] has signaled the beginning 
of extragalactic high-energy neutrino astronomy. After three years of data taking, 36 neutrino events (plus one event 
whose energy and direction cannot be reconstructed) with energies between approximately 30 TeV and 2 PeV have 
provided the evidence for the existence of an extraterrestrial neutrino flux at 5.7 ct [3]. The discovery of this flux 
has motivated a large number of studies in the literature to unravel their origin, from different scenarios within 
standard cosmic-ray sources [4-50] to more exotic possibilities [51-76]. However, the current angular information and 
statistics do not allow the identification of the neutrino sources and the flux remains compatible with an isotropic 
distribution [2, 3]. A larger and improved version of the IceCube neutrino detector has among its major goals the 
resolution of the individual sources of the observed astrophysical neutrinos [77]. However, by focusing on observables 
such as the flavor composition and spectrum, the diagnostic power grows proportionately to the statistics in the 
available sample of events. 

The astrophysical neutrino flavor composition has long been recognized as a powerful tool to disentangle the 
underlying mechanism of ultrahigh-energy neutrino production [78-98], including scenarios of exotic physics [79, 99- 
115]. In the standard scenario, ultrahigh-energy neutrinos are produced by the decay of pions and kaons and secondary 
muons, produced by hadronic interactions in extragalactic cosmic accelerators. This mechanism gives rise to a flux 
flavor ratio at the cosmic source of (ag : : Ofi-js = (1 : 2 : 0)s. Current measurements of the neutrino oscillation 

parameters [116-118] imply that this source composition is transformed into a neutrino flavor ratio at the detector 
position on Earth of {ag : : ai-je — (1:1: 1)0 [119]. Observationally, in the energy range of interest here, 

the neutrino flux flavor ratios are tagged via two possible event topologies: muon tracks, and electromagnetic or 
hadronic showers. Our previous studies of the neutrino flux flavor ratios [120, 121] showed that the canonical ratio 
at Earth, {ag : : Q!r)e = (1:1: 1 ) 0 , assuming the typical E~^ high-energy neutrino spectrum, is disfavored 

at 81% (92%) confidence level (C.L.) by a fit to the two-year (three-year) IceCube data, although when accounting 
for the systematic uncertainties in the backgrounds, the significance is slightly reduced [121]. We pointed out that 
this mild, albeit compelling, tension could be due to several factors, related to a misunderstanding of the background 
events, to a misidentification of track events as showers, or to nonstandard mechanisms of neutrino production at and 
propagation from the sources. 
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In these previous works we only analyzed the total number of events, discarding the additional information encoded 
in both the signal and the background energy distributions. Although the energy behavior has already been exploited 
in the literature as a diagnosis of the origin of the high-energy neutrino events [5, 37, 59, 122-124], it has never been 
used in conjunction with the flavor information and by allowing the background normalizations to vary freely. Spectral 
analyses of the full data have been performed, either fixing the flavor composition to (1 : 1 : 1)0 [2, 3, 124, 125] (see 
also Ref. (37] for other cases) or varying the flavor composition but with fixed backgrounds [123]. In this study, we 
present the first energy and flavor-dependent analysis of the 36 high-energy IceCube showers and track events where 
the astrophysical flavor composition, spectral index and normalization (assuming a power-law flux equal for neutrinos 
and antineutrinos), along with the number of background events are left free to vary. Furthermore, we extend the fit 
performed in our previous works [120, 121] to higher energies, in order to extract crucial information concerning the 
astrophysical high-energy neutrino flux and its origin. 

This study includes the spectral index and the normalization of the incoming neutrino flux, the number of conven¬ 
tional atmospheric muon and neutrino backgrounds (considered separately), and a prompt neutrino flux component 
from charmed particle decays in the atmosphere. We also consider a case in which the extraterrestrial high-energy 
neutrino flux follows a broken power law with a break of one unit in the spectral index at Ei, = 1 PeV. We finally 
examine the effect of a nonzero probability of misidentifying muon tracks as shower events. 

We present the results from our fits to the 36 events measured by the IceCube experiment considering four possible 
energy ranges: the hducial 28 TeV—3 PeV energy range, which covers the publicly available IceCube high-energy 
neutrino events; the 60 TeV—3 PeV range, which was used by IceCube for their inference of the properties of the 
astrophysical flux; and an extension of these ranges to 10 PeV, where a large signal from the well-known Glashow 
resonance [126] is expected — but not seen so far. This serves to illustrate the different observational and physical 
effects that may bias the reconstructed properties of the high-energy neutrino flux. Namely, a larger value of the 
minimum deposited energy would eliminate most of the expected atmospheric muon background events, while a 
larger value of the maximum deposited energy would encompass the Glashow resonance at En ~ 6.3 PeV, which 
should give rise to yet-unobserved events in the few PeV region (see also Ref. [122]). The absence of these events 
could point, among other possibilities, to a break in the high-energy neutrino spectrum around a few PeV, which 
could nicely be connected with the behavior of the ultrahigh-energy cosmic ray spectrum [5, 37, 59, 122-124]. 

Concerning the fiducial 28 TeV—3 PeV energy range, the canonical ratio of neutrino flavor fluxes at the Earth (1:1: 
1)0 becomes more favored once the information from the event energy distribution is included in our multidimensional 
fits. For this case, as well as for the other possible energy ranges explored here, the precise values of the best-fit 
parameters for the neutrino flavor ratios depend on the different quantities involved in the fits. As a general trend, 
the choice of a larger low-energy threshold (from 28 TeV to 60 TeV) implies a wider allowed region in the flavor 
parameter spaces that encompass a signihcant muon neutrino component, even if the best fit shows a slight preference 
for an electron neutrino component. Increasing the upper limit on deposited energy (from 3 PeV to 10 PeV) results 
in a favored tau neutrino component and a steeper spectrum. However, in the latter case, the existence of a break in 
the high-energy neutrino spectrum could restore the preference for a large electron neutrino component. Finally, we 
remark that to firmly establish the neutrino flavor ratios from data and therefore unravel the astrophysical neutrino 
production mechanism, the fraction of tracks that could be misidentified as showers is a critical parameter. Indeed, 
a large fraction of tracks being misidentified would skew the reconstructed flavor ratio away from a significant muon 
neutrino contribution. 

The structure of the paper is as follows. We start in Sec. II by describing our calculations of the shower and track 
neutrino event rates as a function of the deposited energy, including a description of the main ingredients: the atten¬ 
uation and regeneration effects in the neutrino propagation through the Earth; the definition of the electromagnetic 
(EM)- equivalent deposited energy for each type of event in neutrino-nucleon and neutrino-electron interactions and 
how it depends on the geometry of the detector; and our computation of the IceCube effective mass in terms of the 
measured deposited energy. Section III contains a detailed description of the background treatment followed in our 
numerical analyses, whose methodology is presented in Sec. IV. The results arising from our fits to the three-year 
IceCube neutrino data are shown and discussed in Sec. V. Finally, we highlight our most important findings and 
discuss their implications in Sec. VI. 


II. EVENT RATES 

The incoming neutrino fluxes at the IceCube detector are characterized by their energy E^, and their flavor compo¬ 
sition, (cTe : : ar)e) while individual events are characterized via the total EM-equivalent deposited energy Edep 

and the event topology (tracks or showers). Below PeV energies, interactions with nucleons are the dominant ones. 
Showers are induced by both z/g (and i>e and v^) charge current (CC) interactions, as well as by neutral current 

(NC) interactions of neutrinos of all three flavors. Tracks occur when a muon is produced: either via a interacting 
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via W exchange (CC) or via tau lepton decay with the tan produced in a CC event. For each channel, represented 
by a subscript c, the measured deposited energy rate dN"^/dEdep,i, which is evaluated for each observed event i, is 
given by 


dE(\Qpi 



dN^ 

dE^i-^Q 


R{Etrue :-^dep,ij ^(-^true)) ^-^true • 


( 1 ) 


The function i?(i?true, ^'dep.i, ^•(E'true)) Is included to account for the energy resolution cr(£'true) and a fit is provided 
in Appendix B. 

The spectrum of true deposited energies in the detector depends on the attenuation and regeneration factor 
Attl^{E^) due to the absorption of neutrinos when traversing the Earth, the detector effective mass Meff(i?true) 
as a function of the true deposited energy, the incoming neutrino flux d(j)l^{E^)/dE^ of type / (either of astrophysical 
or atmospheric origin) and the production cross sections da^^/dEtme of neutrino Vi for process c. For electron and 
muon neutrinos interacting with nucleons, this can be written in general form as 


dN^ 

dE/true 


= TNa 


Meff(Atrue) 


5 Atrue) 

dEy dF/ti'ue 


dEy , 


( 2 ) 


where T = 988 days is the time of data taking and Na = 6.022 x 10“^^ g“^. 

For tau neutrinos below a few PeV, the produced tau lepton decays inside the detector, so one has to take into 
account the energy spectra of its decay products. In Appendix B we provide the detailed expressions for the differential 
rates for all the processes we consider as a function of the measured EM-equivalent deposited energy, Eq. (1). Let us 
note that we also consider interactions of all neutrino flavors with electrons and we also give the full expressions in 
Appendix B. All the relevant cross sections used in this work are provided in Appendix A. 

Below, we describe how we compute the attenuation and regeneration factors of the different neutrino fluxes in their 
passage through the Earth, we explain how to compute the true EM-equivalent deposited energy of hadronic showers 
and muon tracks and finally, we discuss how to obtain the effective mass which encodes the detector 

efficiency, as a function of the true EM-equivalent deposited energy. The reader not interested in these technical 
details can jump to Sec. III. 


A. Attenuation and regeneration of neutrinos in their passage through the Earth 


The rise of the neutrino-nucleon cross section with energy implies that, for energies above a few TeV, the mean 
free path inside the Earth becomes comparable to the distance traveled [127, 128], so the Earth attenuates the flux 
of neutrinos. The transport equation for vg = {ve,Vfj,,Ve,Vfj,} is given by [129-135] 


dX dEy j 


-Na {a 

+Na f 

Jo 


NC^P I XCfp ^4>tii.Ey,X) 

VI \Ev)EOy^ [Ey)) - — - 

' dy - y),y) dcj,l^{Ey/{l - y),X) 

1- y dy dEy 


( 3 ) 


where (j^^{Ey) and cry^{Ey) are the NC and CC neutrino-nucleon cross sections of (see Appendix A), X{9) = 
Jp p{x) dx is the column depth (with L = 2i?0cos0), the Earth radius, and 9 the nadir angle of the direction 
of the neutrinos with respect to the position of the detector. Throughout this work, the density of the Earth p{x) is 
assumed to be given by the STW105 (also known as reference Earth model, REF) model [136, 137]. 

For the energies of interest for cosmic neutrinos, in general, interactions with electrons can be neglected. However, 
for j?e, the resonant production of a IF boson at energies around Ey = M^/nie ~ 6.3 PeV (the so-called Glashow 
resonance [126]) has a cross section larger than the neutrino-nucleon cross sections and must be taken into account^. In 
all our computations not only do we include this resonant cross section but also all the neutrino-electron interactions. 
Therefore, in the transport equation above, Eq. (3), we make the substitutions 

a^^iEy) + a^^{Ey) ^ a^^{Ey) + a^^{Ey) + {Ey) , (4) 


d(T^p{Ey,y) da^^{Ey,y) ^ dal^^^{Ey,y) 

dy dy dy 


( 5 ) 


^ We do not take into account the small effect due to the Doppler broadening of the resonance [138]. 
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where (see Appendix A) 


^l^Ev) = ^i^^eiEu) + o'p^,;i(Aj/) + ap^ ,^{Ey) + cr|^^^(Ay) , 


( 6 ) 

(7) 

( 8 ) 
(9) 


The notation used throughout this work is as follows: is the cross section for a neutrino of flavor £ to interact via 

I = {electron scattering, CC nucleon scattering, NC nucleon scattering} and produce a lepton £'. In this work, all the 
cross sections are defined per nucleon and we assume the Earth to be an isoscalar medium, so the number of electrons 
is half that of nucleons. This relative factor of 2 is already included in the cross sections above (see Appendix A). 

In the case of electron and muon neutrinos, the charged leptons produced in CC interactions in the Earth are 
quickly brought to rest and either are absorbed or decay at rest, and hence do not contribute to the high-energy flux; 
in the case of and Dr, however, the produced tau leptons can decay before being stopped. Therefore, as long as 
the T decay length is shorter than its mean free path in the Earth, which occurs up to several hundred PeV, I'r ’s and 
Vt ’s are not absorbed, but degraded in energy, so the regeneration of the Vr and 9^ fluxes due to the daughter and 
i>T from T decays must be accounted for [132, 134, 139-144]. For the energies of interest the tau decay length is much 
shorter than the interaction length and r energy losses may be neglected [142, 145-149]. The transport equations for 
I'r and Pr (and tau leptons) are given by the coupled equations [132, 134, 135, 142] 


5A dE J 




+Na 


dy d<C(£;/(i_y),j/) d^l^{E/{l-y),X) 


1-y 


{E/mr)Tp{X) Jo 


dy dE 

^ #A-E'/(1 - y),^) 


dy 


dE 


( 10 ) 


J)_ f dct>f{E,X) \ 
dX \ dE ) 


1 dcpf{E,X) 

{E/mr)Tp{X) dE 


+ NA 



dy da^J^{E/{l-y),y) d(j)l^{E/{l - y),X) 
1 — y dy dE 


( 11 ) 


where r and uIt are the tau lepton lifetime at rest and the tau lepton mass, respectively, d(f)l{E, X)/dE is the flux of 
tau leptons at X{9) produced via CC interactions and dn(z)/dz is the distribution of tau neutrinos after tau lepton 
decays with z = E^^/E^, for which we use the parametrization given in Ref. [140]. Analogously to the case of electron 
and muon neutrinos and antineutrinos, we also include interactions with electrons by making the substitutions, for 
both Vr and Prj 


d(Tj^{E^,y) daJ^{E^,y) dal^^^{E^,y) 


dy 


dy 


dy 


and only for but not for 


where (see Appendix A) 


dcr^J^iE„,y) dajJ^{E^,y) daf,^.,{E^,y) 


dy 


dy 


dy 


<(15.) = <,,(A.)+<,.(A,) , 

^ii-y{Ey) = Ujj JE^) . 


( 12 ) 

(13) 


(14) 


(15) 

(16) 


We additionally consider the secondary r'e and (and v^, and v^) fluxes produced after (and Vr) CC interactions 
and the subsequent r decay into leptonic channels [150, 151], although they are only relevant for very hard spectra. 
To do so, one has to couple the transport equations for ui = with those for i/r by adding the following term 

to Eq. (3) [135, 151]: 


gUe,x) 


_1_ , dur^uA^ - y) d(t)l{E/{l - y),X) 

{E/mr)Tp{X) Jo ^ dy dE 


(17) 
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FIG. 1. Whole-sky averaged attenuation and regeneration factors for astrophysical neutrinos traversing the Earth, 
for an isotropic power-law spectrum E ', with the IceCube best fit for spectral index 7 = 2.3 [3]. Neutrinos (solid lines) and 
antineutrinos (dashed lines) are shown separately, and the attenuation and regeneration factors for Ve, and Ut are represented 
by red, green and blue lines, respectively. Note that the curves for Ue and are almost identical to those for and 15^, except 
in the case of he around the Glashow resonance, which is clearly visible. 


where for fully polarized tau leptons, the distribution of V£ for r decays, with z = is given by [152] 

-0.18 (4-12z-hl2z2-4z3) . (18) 

To compute the primary neutrino flux after attenuation and regeneration and the flux of secondary electron and 
muon neutrinos we have followed the approaches of Refs. [130, 132, 135] and alternatively that of Ref. [153] for he 
(used in a different context), because around the Glashow resonance the former approach presents some issues of 
convergence. We refer the reader to those works for the details of these calculations. 

We compute the neutrino and antineutrino fluxes at the detector for all flavors as a function of the nadir angle 9 
and the neutrino energy Ei, and obtain the 47r-averaged attenuation and regeneration factors for isotropic fluxes, 

= ^ ^1 + J Attl^{E,,, 9) dcos9^ , (19) 

where Att^^{Ey,9) = {d(j)^^{E,L)/dE^)/[d(j)^^{E,Q)/dEy) represents the fraction of the initial flux d(l)^^{E^,0)/dE^ 
propagating through the Earth with nadir angle 9 that reaches the detector. The first term in Eq. (19) represents the 
averaging over downgoing neutrinos that do not cross the Earth. 
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As an example, these whole-sky averaged suppression factors are shown in Fig. 1 for an isotropic power-law spectrum 
with spectral index 7 = 2.3. The factors for neutrinos (solid lines) and antineutrinos (dashed lines) are shown 
separately. One can see the effect of regeneration of inside the Earth (blue lines) and notice that the factors for 
and T'e (red lines) are very similar to those for and (green lines). The only large deviation occurs around the 
Glashow resonance in the case of Pe (red dashed line), as can clearly be seen. 


B. Deposited EM-equivalent energies 

Following the notation of the IceCube Collaboration, we define the true deposited energy Etrue as the total true 
EM-equivalent energy deposited after a neutrino interaction with the vertex within the detector. Each channel has 
different efficiencies when it comes to producing a measured energy deposition in the IceCube detector, since what is 
ultimately measured is the Cerenkov radiation emitted while the charged particles propagate. 

The light yield of the electromagnetic cascades in ice is approximately proportional to the total track length of 
all the electrons in the cascade, which is used to calculate the number of Cerenkov photons and hence, the energy 
of the cascade. On the other hand, the total track length of hadronic cascades is not linear with the energy of the 
cascade. This is due to the presence of more neutral particles like neutrons, to large losses due to the binding energies 
in hadronic processes and to a higher Cerenkov threshold for hadrons [154, 155]. This always results in a reduction 
of the EM-equivalent energy (or total track length) of a hadronic shower with respect to that of an electromagnetic 
shower. Thus, for hadronic showers the EM-equivalent deposited energy is given by 

Eh = Fh{Ex)Ex, ( 20 ) 

where Ex is the shower energy and Fh{Ex) is a suppression factor, which represents the ratio of the track length of 
a hadronic and an electromagnetic shower of the same energy and is provided in Appendix B. 

At these energies, the lifetime of a muon is much larger than the time it takes to cross the detector. The total 
energy deposited by a muon can be described by the mean stopping power. Although at the energies of interest these 
energy losses can be stochastic and large fluctuations around the mean are expected, it is reasonable to treat them as 
continuous and approximate the muon energy loss rate by 

(-'^'^=a{E^) + h{E^)E^ (21) 

where E^ is the muon energy, a{E^) is the electronic stopping power and b{E^) takes into account radiative processes 
(bremsstrahlung, pair production and photonuclear interactions). Writing the mean stopping power in this way is 
convenient, since both a{E^) and b(E^) vary slowly in the energy range of interest here. For the sake of simplicity, 
in this work we use a fit obtained from tabulated data for the muon loss rate in ice [156], which is provided in 
Appendix B. 

To compute the energy deposited along a given muon track, we need to know the position where the muon was 
produced and its direction, i.e., the interaction vertex in the detector as well as the distance from the vertex to the 
edge of the detector volume. Bearing in mind that this deposited energy is, in general, much smaller than the energy 
deposited by the hadronic shower, we compute the vertex position and direction-averaged energy deposition, (AE^), 
along an average muon track in the detector given the initial muon energy, To do this, we approximately describe 
the detector volume as two stacked cylinders of radius 500 m, with respective heights 275 m and 545 m, separated by 
an 80 m dust zone [2]. We consider all the points and directions in the detector to be equally likely and do not count 
in any muon produced in the dust zone, although we do take into account energy losses in that zone. We illustrate 
this geometry in Fig. 2. 

For a given interaction vertex location in cylindrical coordinates (r, z) inside the total detector volume, the angle- 
averaged distance traveled by a muon inside the detector is 

1 p2-R / ^cos 6»o(i-,z,0) \ 

{L{r,z)) = — / d(/) / £rir, z, 6, (l))d cos 6 + / 4(r, z, 0, ())) dcos 0 , (22) 

Jo yJ0 Jcos 9Q{r,z,(p) J 

where Oq is the angle to the intersection between the floor and wall of the cylinder and ir is the propagation distance 
from the vertex to the wall if 0 > 9q, and 4 is the distance from the vertex to the floor, when 9 < 6q. Upgoing 
and downgoing tracks produce an equal contribution to {L{r,z)). Therefore, the integrals in Eq. (22) cover only 
downgoing muons, while an overall factor of two ensures that the full angular range is considered. The full average 
over all interaction vertex points is then 


(L) = :^ y dz J drr {L{r, z)) , 


(23) 



FIG. 2. Approximate geometry of the IceCube detector, which we use to compute the average distance traveled by a muon 
generated inside the detector volume in Eqs. (22)-(23) and thus the expected energy deposition of a muon track, Eq. (25). From 
an interaction vertex located at {r,z), the distance £ represents the distance traveled by a muon propagating in the direction 
(6^, (j>). When 9 < 6q, £ = £z and the muon exits via the detector floor; when 9 > 9o, £ = £r and the muon exits via the “wall”. 
We ignore events generated in the dust zone (gray band), but we do include the energy lost in that zone. 


where the volume Vtot covers each point in the cylinders, outside the 80 m dust zone. Numerically integrating Eq. (23) 
yields an average distance of (L) = 406 m. 

From Eq. (23) it is easy to compute the average muon energy loss per track, (AE^), but substituting the lengths 
£r (or ^z) by 

e{r,z} ^ 1 - , (24) 

where b is defined in Eq. (21). Evaluating this integral yields the total averaged deposited energy along an average 
muon track, in terms of the initial muon energy E^, 

{AE^} = (1 - [E^ + a/b) = E^ (E^ + a/b) , (25) 

where the factor for the average fraction of energy lost along a track of a muon produced at the neutrino interaction 
vertex is found to be = 0.119. 

Finally, the detector geometry is also important in the case of r lepton production. When the t energy exceeds a 
few hundreds of TeV, its lifetime is long enough that a significant fraction of tau leptons would escape the detector 
volume before decaying. In a similar way as done for muons, the average fraction of taus with energy E^. decaying 
inside the detector is given by Eq. (23), but with the substitution 

(■{r,z} 1 - . (26) 

Here we are neglecting the tau lepton energy losses, which scale inversely proportional with the mass and thus are an 
order of magnitude smaller than for muons [145] . Therefore, we can write the fraction of tau leptons decaying inside 
the detector as 


DriEr) = (1 - . 


(27) 


Once the tau decays inside the detector, it has an ~ 18% probability of producing a muon. The average energy 
loss of such a muon is given by the average energy loss, Eq. (25), weighted by the fraction of taus that decay inside 
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the detector, 

= ((1 - (l - ) (E^ + a/b) = {E^ + a/b) , (28) 

where the factor for the average fraction of energy lost along a track of a muon produced in a tau decay after a 
Vr or Dr CC interaction inside the detector is given by Ef^^riEr)- Fits to both, Dr{Er) and E^^r{Er), are given in 
Appendix B. 


C. IceCube effective mass 


We now turn to the computation of the effective mass Meff(Atrue) in terms of the deposited energy in the detector, 
rather than the neutrino energy, since the former is the actual observable. This effective mass can be seen as the mass 
of the target material times the efficiency of converting (EM-equivalent) deposited energy into an observed signal. 
Obtaining this is not completely straightforward, since the effective masses provided by the IceCube Collaboration 
are given as a function of the incoming neutrino energy E^,, for neutrino-nucleon NC and CC interactions involving 
the three flavors. To perform analyses in terms of the differential deposited energy spectrum, as we present below, we 
need to express the effective mass in terms of the quantity Atrue and thus, we need to perform a deconvolution. Using 
the fact that the effective masses Ma{E^) provided by the IceCube Collaboration are meant to be postinteraction, 
they are related to the effective mass Meff(Etrue) via 


««(f;.) + <«(a.)) 

{a^^{E^)+a^^{E^)) 


i dy +—y - ) 

^'fCC(F) (d<jZ^{K,y) dag^{E,,y)\ 

= Jo [ dy + dy ) 



d^ug {Ey,y) da^g iE^,y)' 


dy 


dy 


M,s{E^^)dy , 

(29) 

M,s{E^^)dy , 

(30) 

M,BiEY)dy , 

(31) 


[cr 


.CC 


{Eu) +1 


Jo \ dy dy J 


xf Y. {Dr{Er)M,siEg;^) + {l 

® k—h,e,fj. 


Dr{Er)) M,e{Eh)) 


'^d,dy{32) 

dz 


where dnk{z)/dz [140] is the energy distribution of the daughter Vr, e, or p with energy E^^^r (z = E,y^^r/Er), E^^r 
{z = Ee^r/Er), Or [z = E^^r/Er) from T decay via the hadronic, electronic, or muonic channel {k = {h, e, ^}), 
respectively, and the true EM-equivalent deposited energies for each case are given in Appendix B. Note that the 
function is common to all flavors and channels, since it only depends on the true EM-equivalent deposited 

energy Etrue- However, notice that Etrue is not given by the same expression for the different types of interactions. 
Let us also point out that we do not include the energy deposition along tau tracks, which is negligible, as the losses 
in this case are much smaller than in the case of muon tracks. 

To obtain Mefj(i?true)i we need to find a suitable parametrization for it. We consider a simple functional form with 
three independent parameters, given by 


if x > 0 


Meff(x) = 



where pice = 0.9167 g/cm^ is the ice density and 


X = log 


10 


Ef 


Ath 


with Eth = 10 TeV. 

Then we perform a least squares analysis for the three parameters {c, d, q}, 


(33) 


(34) 


X 


2 


E 


MfgE,,,;c,d,q) \ 


2 


(35) 
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FIG. 3. Effective masses as a function of the neutrino energy. Solid lines are the effective masses provided by the 
IceCube Collaboration, whereas dashed lines represent the effective masses computed with Eqs. (29)-(32) using the best-fit 
values for {p, q, c} for the effective mass as a function of the deposited energy, Eq. (33). Red, green, blue and gray histograms 
represent CC Ue, Ufj,, Vr and NC (all flavors) interactions. The black histogram depicts the effective mass function Mefl(i?true), 
Eq. (33), as a function of the true deposited energy, i?true. 


where i indicates the neutrino energy bin (we use all bins where the effective mass is different from zero) and 
j = {NC, — CC, — CC, — CC} indicates the type of interaction. The actual effective masses are 

taken from Ref. [2] and Mj^{Ej,y,a,b,c) are obtained from Eqs. (29) - (32). Up to two significant digits the best-fit 
parameters are found to be c = 0.50, d = 1.1, and q = 4.6. 

The accuracy of the fit can be seen in Fig. 3, where we show the IceCube effective masses with solid lines and our 
reconstructed results, computed with Eqs. (29)-(32) using the best-fit values for (c, d, q} for Meff(Etrue), Eq. (33), are 
shown with dashed lines. The result for CC CC, NC) interactions is depicted in red (green, blue, gray). 
We also illustrate the effective mass as a function of the deposited energy. As can be seen, the accuracy of the fit is 
very good over the whole energy range, and, in particular, below the 10% level for the relevant energy range. We note 
that the resulting mass as a function of the deposited energy is very similar to the effective mass for (plus i>e) CC 
interactions as a function of the neutrino energy. This is expected, as the energy deposited after Vf. CC interactions is 
very close to the actual neutrino energy (one electromagnetic shower that fully deposits its energy and one hadronic 
shower that deposits most of its energy). We also point out that the effective mass for CC interactions is larger 
than for CC below a few hundred TeV, due to the fact that the produced tau lepton in the former case, after 
decaying, deposits more energy than the produced muon in the latter case. There is also an interesting feature when 
comparing these two histograms. At energies close to a PeV, they cross each other. This has to do with the fact that 
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at these high energies some tan leptons exit the detector, which decreases the efficiency for this type of interaction. 
Finally, the effective mass for NC interactions is the smallest one because only one hadronic shower could be detected 
and an important fraction of the neutrino energy escapes the detector in the form of daughter neutrinos. 


III. BACKGROUNDS 


The sources of background for cosmic neutrinos in IceCube are due to the interactions of cosmic rays with the 
nuclei of the atmosphere. These interactions produce fluxes of secondary muons and neutrinos from all directions, 
the so-called atmospheric muon and neutrino fluxes. 

Although the rate of atmospheric muons is huge, about 3 kHz in IceCube, it can be significantly reduced by using 
part of the detector as a veto for entering muon tracks [2]. The final event selection allows a reduction at the level of 
one part in 10 ® for muons with at least 6000 photoelectrons detected (approximately, deposited energies higher than 
30 TeV), so the final expected event rate is of a few events per year. This is estimated based on data by tagging 
muons in one region of the detector and by measuring the detection rate in a separate region equivalent to the veto 
region. The expected number of events above the threshold after 988 days of data taking is 8.4±4.2. Nevertheless, to 
perform a statistical analysis based on the energy distribution of the detected events, we need to know the spectrum 
of the veto-passing atmospheric muons. To good approximation, the convolution of the atmospheric muon spectrum 
and the veto probability is well described by a power law in the range of study [157]. To determine the spectral index 
of this distribution we use the IceCube expected rates in two energy bins [3]. For deposited energies below 60 TeV 
the expected number of muons after 988 days is 8.0, whereas for energies above 60 TeV, only 0.4 atmospheric muons 
are expected to contribute to the background. Therefore, the spectral index 7 ^ of the distribution is given by^ 


^ log(2I) 
log(60TeV/F;„i„) 


(36) 


where is the minimum measured EM-equivalent deposited energy considered in the sample. However, the 

actual threshold is given in terms of the minimum number of detected photoelectrons (6000), not in terms of the 
minimum deposited energy. The conversion from the number of detected photoelectrons to deposited energy is not 
straightforward, since the scaling factor depends on the region of the detector where each event occurs. However, 
a detailed and careful analysis in terms of detected photoelectrons requires the full Monte Carlo that describes the 
detector and, therefore, it has to be performed by the IceCube Collaboration. Here, we use Emin = 28 TeV, which 
is the round value below the energy of the less energetic of the events in the full sample. Let us note that we have 
checked that small variations in this minimum energy do not result in significant changes in our final results. For 
this value of Emim 7 /i — 5. On the other hand, the misidentification of tracks as showers for the muon atmospheric 
background is quoted to be < 10% [3]; we take it to be exactly 10%. This choice does not have a significant impact 
on the results. 

In addition to the veto-passing atmospheric muons, atmospheric neutrinos constitute the other source of background 
in this search. In this work, we consider the conventional (mainly from tt and K decays) and prompt (from charmed 
meson decays) atmospheric neutrino flux calculations of Refs. [158-161] based on the Hillas and Gaisser (HGm) cosmic- 
ray approximation [162, 163] and the hadronic model of Kimel and Mokhov (KM) [164], with updated parameters [165- 
167], for the conventional flux, and on the Zatsepin and Sokolskaya cosmic-ray model (ZS) [168] and the quark-gluon 
string model (QGSM) [169-172] for the prompt flux. We use the (conventional and prompt) atmospheric i/g, 

Vf, and fluxes as a function of the nadir angle and computed up to 100 PeV. This calculation results in similar 
spectra to those computed by Ref. [173] in the case of the conventional flux, which the IceCube Collaboration uses. 
However, the latter calculation only extends up to 10 TeV, so it has to be extrapolated up to PeV energies [174]. Let 
us note that the ZS±QGSM prompt flux results in a larger number of expected events than the flux obtained with 
the dipole model [175, 176], although important differences in shape are only present above PeV energies due to the 
ZS cosmic-ray model. Nevertheless, at those energies prompt neutrinos are not expected to have any impact on the 
discussion. 


Although it has no impact on our statistical analyses, to plot event spectra in Sec. V, we assume a broken power law with breaking 
energy at 28 TeV, such that the number of atmospheric muons in the range 10^’^ TeV—10^ ‘^ TeV after 988 days in IceCube is 0.2 [3]. 
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In this analysis we also include the suppression of the flux of downgoing atmospheric t'^’s and P^’s {f = v for 
conventional and f = p for prompt^) contributing to the background by tagging the muon produced by the same 
parent meson decaying in the atmosphere which can trigger the muon veto of the detector. To compute the passing 
probability we have followed Refs. [177, 178]. This probability is defined as [178] 


pUE,,e,) = i- 






6^)/dE^,d cos 9^ 


(37) 


where d(j)l’*{E^,6^,Efj,^^in)/dE^ d cos 9z is the flux of neutrinos that reach the detector with zenith angle 9z = it — 9 
and are accompanied by the muon produced by the decay of the same meson, and it is given by 


where 




dEi, d cos 9z 


1 — Z]^N 


i—TT,K 


Bvi 


1 + Bi El, cos 9*/e 


Ai = 


Zn,i 1 


1 - 1 + 7i- 


,7.^ + 1 ’ 


Bi = z,- 


7 i/ + 2 1 — Ajv / Ai 


7i, + 1 IuAj/Aat 


1 F, 

= max{--,1 + —} , 


1 — r. 


E„ 


m: 


2 ’ 


(38) 


(39) 

(40) 


with 7 ^ = 1.7, Znn = 0.298, An = 120 g/cm^, A^ = 160 g/cm^, Ak = 180 g/cm^, e,r = 115 GeV, ck = 850 GeV, 
Br^ = 1, Br^ = 0.6355, Zn^-k = 0.079 and Zn^k = 0.0118 [163, 179]. The minimum muon energy at production that 
is required for the muon to reach depth Xjc = die/cosmic with at least energy Eth = 10 TeV is 


E^i.^inin = Et,i " l) ^ ’ 


(41) 


where a and b are given in Eq. (21), we take die = 1.45 km as the vertical depth of IceGube and 


cos die 



sin^ 9z 

(1 + die/ {R(b - dicW 


(42) 


is the local angle at the surface of the Earth with respect to the position of the detector taking into account the 
curvature of the Earth. 

The quantity cosd* is the cosine of the local zenith angle at the altitude of the first interaction^, datm = 32 km, and 
takes into account the curvature of the Earth, representing a non-negligible correction for angles larger than 70° [152]. 
Analogously to Eq. (42), 


cos 



sin^ 9z 

(l+/iatm/l?®)" ■ 


(43) 


Finally, to compute Eq. (37), the parametrization for d(j)l^(E^,9z)/dEydcos9z is given by Eqs. (38)-(43), but with 
Zi,min = 1/(1 — ri). Moreover, as done by the IceGube Collaboration [2, 3], the suppression factor, i.e., the second 
term in the right-hand side in Eq. (37), is bounded from above at 0.9 to cover uncertainties in hadronic interaction 
models. 

Once we have the passing probability for atmospheric and (conventional and prompt), we can compute the 
whole-sky averaged attenuation and regeneration factor, which is given by 


{Ef ) — 


CPiiE.,0) 


AE.,e) 


dEiyd COS 9 


d cos 9 + /p Attl {E,z,9) 


dEiydcos 9 


dcosd 


/- 


dE„,e) 


1 dEjyd COS 9 


dcosd 


(44) 


^ Rigorously, one should compute the veto-passing probability for prompt atmospheric neutrinos using their corresponding fluxes. We 
simply use the same result as that for conventional neutrinos. 

Although the average altitude of the first interaction depends on Oz, we take its value at Oz ~ 7r/2 [180]. 
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FIG. 4. Whole-sky averaged attenuation and regeneration factors for atmospheric neutrinos traversing the Earth, 
Attl^iE,,), using the HGm+KM model (red lines for i/e and Ue and dark green lines for and Vff) for the conventional neutrinos 
(/ = u) and the ZS+QGSM model (orange lines for and i>e and light green lines for and Vff) for prompt neutrinos from 
charmed meson decays (/ = p), from the calculations of Refs. [158-161]. Neutrinos (solid lines) and antineutrinos (dashed 
lines) are shown separately. 


On the other hand, for the nadir (or zenith) angle-dependent (conventional and prompt) atmospheric and fluxes, 
the corresponding 47r-averaged attenuation and regeneration factor is 


Attl (E^) = 


pO d4>l {Ely,9) , /I , Axj-f (f\\ j n 

f-1 dE^d cos 9 + /o dE^d cos 9 ^OS ^ 


/- 


d(pl^{Eiy,9) 
1 dE,yd cos 9 


dcosO 


(45) 


In both cases, we have explicitly written the nadir (or zenith) angle dependence, which we use in our calculations. 
The results, using the calculations of Refs. [158-161], for the conventional atmospheric neutrinos with the HGm-|-KM 
model (red lines for Vg and Dg and dark green lines for and and for the prompt atmospheric neutrinos with the 
ZS-I-QGSM model (orange lines for t'g and Pg and light green lines for and P^), are shown in Fig. 4. Again, we 
depict the suppression factors for neutrinos (solid lines) and antineutrinos (dashed lines) separately and the Glashow 
resonance is clearly visible for both types of fluxes. Note that the attenuation and regeneration factors for prompt 
and Pg atmospheric neutrinos are similar to the factors for astrophysical Vg and Pg neutrinos shown in Fig. 1, unlike 
what happens for conventional neutrinos and for and P^ prompt and conventional fluxes. The former is due to the 
fact that the conventional flux is nonisotropic, with the maximum around the horizon, where very little absorption 
takes place. The latter is due to the extra veto applied to and P^ atmospheric neutrino fluxes, as discussed above. 
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IV. STATISTICAL ANALYSIS 


We perforin an unbinned extended maximum likelihood analysis, using the energy and the event topology informa¬ 
tion of all of the 36 high-energy IceCube events. The full likelihood is given by 


L = 




m 




n n 


2=1 


2=1 


(46) 


where ^v, and Np refer to the expected astrophysical (signal) neutrino, conventional atmospheric neutrino, 
atmospheric muon and (when used) prompt atmospheric neutrino total number of events, respectively, and and 
-^obs total number of observed showers and tracks, all in the energy interval under analysis. The partial 

likelihood for each type of event topology, k = {tr, sh}, is 

= Na iPf({a}, 7) + + Np + Np , (47) 

where for each event observed at IceCube, we obtain the probability density function (PDF) for each type of neutrino 
flux and dehne the normalized PDF, Vf’^, as the probability distribution for an observed event with energy ifdep,* 
and topology k caused by a flux of type / of incoming neutrinos. For an isotropic flux of astrophysical neutrinos 
with flavor combination at Earth {ug : : Or}® and spectrum dcff"/dE^ oc (assuming the same spectrum and 

normalization for antineutrinos), 


n''’“(W,7) = 




E„ 

Ern 


dE, 


deP d£de: 


E 


dNf 


ag 


dE, 


dep ,2 


(48) 


The integral in the denominator covers the EM-equivalent deposited energy region considered for each analysis and 
the sum goes over the three neutrino flavors, i = {e, /r, r}, and the type of event topology, j = {tr, sh}. The spectra 
dN^^’°‘ /dEdep and dN^'^’°‘/dE^ep resulting from the sum of all the partial contributions from different processes to 
showers and tracks from neutrinos and antineutrinos of flavor I are detailed in Appendix B. 

Likewise, for the atmospheric backgrounds we have 


Vi 


' i 


'pk,P 
' i 


1 dN^'' 

1 

Y' l-En,^^ dNi’!^ dEdep i 

^3 JE^,^ '^■^dep d®P.* 

1 dN^’P 


E r^'max 7771 dN^ 


(49) 

(50) 

(51) 


where the conventional and prompt atmospheric neutrino fluxes are used to compute the event distributions just as in 
the astrophysical case. However, note that for the conventional and prompt atmospheric neutrino fluxes the relative 
contributions from each flavor, and for neutrinos and antineutrinos, are obtained by using the corresponding flux for 
each case. The event distribution of atmospheric muons was discussed in Sec. III. For the three background sources, 
we write the sum over flavors of the (neutrino-plus-antineutrino-induced) differential event spectra. 

We perform our analyses of the IceCube events using the MultiNest nested sampling algorithm [181-183], with 
the likelihoods given above. We let the following parameters vary freely: 


{ae,ap,^,Na,N„,Np} , (52) 

where ar = 1 — ag — and Np = 0. We refer to this parameter set as “6P”. We compare these results with 
MultiNest runs where we hx Ni, and Np to the expected values provided by IceCube: = 6.6 (2.4 when the 

energy range is restricted to Edep > 60 TeV) and Np = 8.4 (0.4 above 60 TeV), Np = 0, and the energy power-law 
index of the astrophysical flux to their IceCube best-ht value, 7 = 2.3 [3]. This parameter set is called “3P”. 

To probe the effect of a possible prompt atmospheric contribution from the decay of charmed mesons, we define 


7P = {ag,ap,^,Na,N^,Np,Np} , (53) 

with the further imposition of a half-Gaussian prior on the number of prompt atmospheric neutrinos Np widths 
o'Np = 5.5 (3.2) for the > 28 TeV (> 60 TeV) energy ranges. This reflects the Icr C.L. exclusions provided by the 
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lower energy fluxes quoted in Ref. [3]. For consistency, this case also includes a continuous Poisson prior on 
centered at the expected fluxes given above. 

Finally, for the cases of a broken power-law spectrum and when including track misidentification, we also show 
results varying and setting the other parameters to their IceCube best-fit values. These analyses are 

referred to as “4P”. 

For our analyses we allow the full range of flavor compositions as seen at Earth: may go from 0 to 1 for each 

flavor £ = {e,/r, r}, as long as their sum is unity. In reality, neutrino oscillations are averaged out after propagation 
over large distances, leading to a much smaller allowed region in flavor space as the flavor composition averages to a 
small part of the parameter space at Earth, around (1 : 1 : 1 ) 0 . If the origin of the high-energy neutrinos is purely 
charged pion decay, the expected composition (1:2: 0)5 would oscillate to an observed ~ (1.04 : 0.99 : 0 . 97)0 ratio, 
while neutron decay sources yielding (1:0: 0)s give rise to a ~ (0.55 : 0.24 : 0 . 21)0 flavor ratio at Earth. For 
reference, the full space of allowed source neutrinos is shown in our ternary plots as a thin blue triangle. At present, 
the paucity of observed events gives a very flat likelihood within this area, resulting in very little statistical power with 
regards to the IceCube neutrino composition at the sources. Furthermore, we find that the Bayes factor comparing 
the full flavor range to the averaged-oscillation triangle is close to one, meaning that there is no evidence at present 
for such a restriction to be particularly good or bad. 

We compare the results with different choices for the minimum and maximum EM-equivalent deposited energies 
for the analysis, Emin and Emax, and show how these choices affect what one can conclude about the data sample. 
We consider four distinct energy ranges: 28 TeV—3 PeV, the range which covers all published high-energy IceCube 
events, 60 TeV—3 PeV which eliminates most of the background atmospheric muon events and is the energy interval 
considered by the IceCube Collaboration to present their results, as well as 28 TeV—10 PeV and 60 TeV—10 PeV. By 
extending the analysis to 10 PeV, we cover the energy region of the Glashow resonance at E^, ~ 6.3 PeV where a few 
events are expected if the astrophysical spectrum extends beyond a few PeV. The absence of any observed events above 
2 PeV could indicate a small electron neutrino component, a very soft spectrum, or a break in the spectrum around a 
few PeV. This could nicely be connected with the ultrahigh-energy cosmic-ray paradigm [13, 21, 37, 41]. Conversely, 
the main effect of removing events below 60 TeV is a vast reduction in the atmospheric event contamination. While 
this yields a much cleaner astrophysical signal, it also removes 16 of the 36 observed events, reducing the statistical 
power of the analysis. 

To avoid subtleties related to the problem’s particular geometry in the large parameter space, we present most 
of our results as profile likelihoods of the plotted parameter space, rather than as Bayesian posteriors. To do this, 
we bin the samples produced by MultiNest with respect to the relevant parameters, and within each bin j in the 
parameter space, we find the point of maximum likelihood. Then, we define the \j = —2 log(£j/£niax), where £max is 
the likelihood of the overall best-fit point. The test statistic is assumed to be distributed like a with n degrees 
of freedom, with n equal to the number of parameters in the corresponding analysis. When instructive we also present 
quantities obtained with the Bayesian posteriors. 


V. RESULTS 


In Table I we show the best-fit points from our analyses, as well as the Bayesian posterior means. In the following 
subsections, we present our results and summarize the conclusions that may be drawn from our analyses of the 
high-energy IceCube events. In Fig. 5 we show the spectra of the best-fit points for the 6 P analyses, divided into 
shower-only (top), track-only (middle), and total (bottom) contributions, along with the binned IceCube data with 
Feldman-Cousins errors [184]. The bottom panels also show the atmospheric neutrino and muon contributions to our 
best fit for the cases with an upper cut at 3 PeV. On the left we show the spectra of the different contributions for 
events observed above 28 TeV, while the panels on the right show fits to events above 60 TeV, where the background 
rates are much lower. On the right panels we also show the result for the IceCube best fit in the energy interval 
60 TeV—3 PeV, which assumed (1:1: 1)0 and obtained 7 = 2.3 and N^f, = 4.5 x 10“^® GeV“^ cm“^ s“^ sr“^, where 
is the total astrophysical flux at 100 TeV. From these plots, one can see how the events below 60 TeV push the 
astrophysical flux toward a softer spectrum. Let us also stress the importance of the flavor composition, by noting 
that cases with similar total spectra in the range where they are fitted, give rise to different shower and track spectra 
that allow the breaking of the degeneracy. 
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FIG. 5. Spectra of the best fits of our six-parameter (6P) analyses, overlaid with the best fits provided by the IceCube 
Collaboration and the binned high-energy neutrino event data with Feldman-Cousins errors [184], Left panels: the spectrum of 
the 36 events above 28 TeV. Right panels: the spectrum of the 20 events detected above 60 TeV. Top panels show the shower 
component only; middle panels show the track component, and bottom panels show the total number of events (showers -|- 
tracks), in addition to explicitly showing the atmospheric event contributions to the 28 TeV—3 PeV (60 TeV—3 PeV) case in 
the left (right) panels. 
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FIG. 6 . Left panel: contours in the 7 - plane, at Icr C.L. (dark colors) and 2(7 C.L. (light colors). Filled contours 
(closed curves) represent the cases with a high-energy cut in the deposited energy of 3 PeV (10 PeV). Including the full 
energy range yields a preference for a larger spectral index due to the lack of events from the Glashow resonance. The blue 
star (diamond) is the best-fit point, = 2.6 x 10“^® GeV“^ cm“^ s“^ sr“^ {N^ = 5.6 x 10~^® GeV“^ cm“^ s“^ sr“^), for 
60 TeV—3 PeV (60 TeV—10 PeV), while the red star (diamond) is the best-fit point, = 6.9 x 10“^® GeV~^ cm“^ s“^ sr“^ 
{N,j, = 7.7 X 10“^® GeV“^ cm“^ s“^ sr“^), for 28 TeV—3 PeV (28 TeV—10 PeV). The black square is the IceCube Gollaboration’s 
quoted best fit, = 4.5 x 10“^® GeV“^ cm“^ s“^ sr“^ and 7 = 2.3, whereas the black circle is their best fit when the spectral 
index is fixed to 7 = 2, = 2.85 x 10“'^® GeV“^ cm“^ s“^ sr“^ [3]. Right panel: p values along the range of astrophysical 

events Na. The results in both panels refer to 6 P analyses. 


A. Spectral index 


By letting the power-law index 7 of the astrophysical neutrino flux vary freely, for the 6P analysis we find a best 
fit 7 = 2.96lg'37 for the energy range 28 TeV—3 PeV. When events below 60 TeV are omitted, this is lowered to 
7 = 2.34^0 3®, in perfect agreement with the best-fit value that is obtained by the IceCube Collaboration [3]. This is 
an indication of contamination by the steeper background spectra (subdominant above 60 TeV), which can explain a 
large fraction of the events below 60 TeV and force the astrophysical spectrum to be softer. However, this results in 
a lower probability for events above a few hundred TeV. Conversely, preference for a steeper spectrum when energies 
above 3 PeV are included can be interpreted as an indication of the missing events due to the Glashow resonance. 
We discuss this more thoroughly in Sec. V C. Let us also note that a different analysis, with a low-energy threshold at 
1 TeV, also obtained a steep spectrum with 7 = 2.46 as its best fit [185]. A potential problem of such a soft spectrum, 
for hadronuclear scenarios where neutrinos are produced from pion decays along with gamma rays, is the violation 
of the isotropic gamma-ray background bounds when extrapolated to lower energies [12]. Such a limit would imply a 
power-law index 7 < 2.2 and could in turn help distinguishing different mechanisms of neutrino production. 

In the left panel of Fig. 6 we show the Icr and 2a contours (for 2 degrees of freedom) in the 7 - plane, where 
the total flux normalization in units of 10“^® GeV“^ cm“^ s“^ sr“^, is defined as 


^ = nJ 

dE^ \ 100 TeV 


(54) 


We show with filled contours (closed dashed curves) the two cases with Amax = 3 PeV (10 PeV). The contours for the 
low-energy cut at Amin = 28 TeV (60 TeV) are depicted in orange (purple). We also show the best-fit combinations 
of ( 7 , N,p) inside our contours, indicated by stars (diamonds) when the highest deposited energy in the analysis is 
3 PeV (10 PeV), as well as the best-fit values quoted by the IceCube Collaboration for fixed spectral index 7 = 2 
(black circle) and for their best fit, 7 = 2.3 (black square) [3]. The IceCube analysis assumed the flavor combination 
expected from pion sources, i.e., (1:1: 1)0, and considered the energy range 60 TeV—3 PeV, so it has to be compared 
with our purple contours. The major difference comes from the fact that we do not fix the flavor composition, but 
let it freely float. It turns out that both IceCube points are within the la C.L. contour, although close to the edges. 
From this figure, one can see the behavior, just pointed out, that extending the energy range beyond the window 
considered in the IceCube analysis has a significant impact on the inferred spectral index of the astrophysical neutrino 
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FIG. 7. Effect of accounting for systematic errors. Ternary plots of the profile likelihood exclusions of the neutrino flavor 
composition of the 36 high-energy events seen at IceCube after 988 days, in the energy range 28 TeV—3 PeV. Left panel: 
varying the full parameter space ( 6 P analysis). Right panel: fixing the spectral index to 7 = 2.3 and the backgrounds to the 
expected total number of events given by the IceCube Collaboration [3]: Nf^ = 8.4, Ni, = 6.6 and Np = 0 (3P analysis). The 
black (cyan) lines represent the Icr {2a) C.L. allowed regions. The best fits, indicated by the white circles, are (0.75 : 0.25 : 0)0 
(left panel) and (0.92 : 0.08 : 0)© (right panel). The thin blue triangle represents the space of allowed astrophysical neutrinos 
assuming averaged oscillations during propagation from the sources and the star represents the canonical ( 1 : 1 : 1 ) 0 . 


flux: adding events below 60 TeV (orange regions) or above 3 PeV (dashed lines) steepens the spectrum. This was 
already noted in an independent analysis of the IceCube data [37]. 

When one adds the potential contribution of atmospheric neutrinos from charmed meson decays [7P analyses; see 
Eq. (53)and the description of priors below it], the best-fit astrophysical flux tends to be slightly harder in all cases. 
This mainly has to do with the fact that the background from atmospheric muons is forced to be smaller than in the 
6P case (where no priors are applied), and hence, a slightly harder astrophysical spectrum is needed to account for 
the number of observed events. 

In the right panel of Fig. 6 we show the p value for the total number of astrophysical neutrinos in the four energy 
intervals. Above 28 TeV we expect ~ 26 — 27 of the 36 observed events to be of astrophysical origin, but with a large 
la spread of ~ 9 — 10 events. Above 60 TeV, the number is ~16 out of 20, with an uncertainty of ^5 events. The 
larger number of background events obtained for the best fit in the 60 TeV—3 PeV case reflects in a slightly smaller 
number of total events from astrophysical neutrinos. Nevertheless, this is not statistically significant at present. 

When we perform the 3P analyses, i.e., fixing the spectral index and the atmospheric background to the IceCube 
analysis best-fit values, our best-fit number of astrophysical events for the interval 28 TeV—3 PeV is reduced to 
~ 20 — 21. This is due to the fact that by fixing the spectral index to 2.3 a larger number of background events are 
needed to explain the data. Otherwise, this can be compensated by a steeper spectrum, as can be seen from Table I. 
For the cases with a low-energy cut at 60 TeV, there is basically no difference because the value of the fixed spectral 
index is very close to the best fit. 


B. Flavor composition 

As was concluded from the results of our previous analysis in the 28 TeV—3 PeV range [120, 121], which omitted 
spectral information and included some simplifying assumptions, a significant component in the astrophysical flux 
is disfavored (see, however. Sec. VD). This is due to the paucity of tracks in the observed event sample. We show 
here that conclusions about the flavor composition at Earth (oe : : Or)® of the high-energy events are strongly 

dependent on the energy interval considered for the analysis and whether the spectral index of the astrophysical flux 
and the number of background events are allowed to vary freely. 

We first consider the range 28 TeV—3 PeV, which is shown in Fig. 7 for the case of the 6P (left panel) and the 3P 
(right panel) analyses. From the right panel, we note that our previous results [120, 121], using a single energy bin, 
fixing the backgrounds to the IceCube expectations and only using the event topology information, is qualitatively 
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FIG. 8. Effect of extending the energy interval of the analysis up to 10 PeV. Ternary plots of the profile likelihood 
exclusions of the neutrino flavor composition of the 20 high-energy events seen at IceCube above Edep = 60 TeV with the 6P 
analyses, in the energy ranges 60 TeV—3 PeV (left panel) and 60 TeV—10 PeV (right panel). The best fits are (0.98 : 0.02 : 0)0 
(left panel) and (0.01 : 0.01 : 0.98)© (right panel). The consequence of not having observed any event above 2 PeV is the shift of 
the best-fit point from a flux with a dominant Vf. component (for i?max = 3 PeV) to an almost pure Vr flux (for Smax = 10 PeV). 
Same format as Fig. 7. 


confirmed when the energy distribution is included (see also Refs. [123, 124]), although adding the spectral information 
as we do here ameliorates the fit for (1:1: 1 ) 0 , bringing the p value of that point, for 7 = 2.3, from 0.14 to 0.29. 
When the backgrounds and the spectral index (left panel) are left free to vary, the canonical (1 : 1 : 1)0 scenario is 
well within the Icr C.L. contour and strong conclusions cannot be drawn. The expected number of background events 
is allowed to be sufficiently small so that a significant contribution to the track sample from astrophysical neutrinos 
is required. Consequently, this also implies a larger number of the total number of astrophysical events. This is not 
surprising, since the disfavored astrophysical component arose from the large number of tracks expected from the 
atmospheric muon and neutrino backgrounds, ^12 in the 28 TeV—3 PeV deposited energy range, compared with the 
actual 8 tracks observed after three years. At the same time, as discussed above, a steeper spectrum is preferred. 

Moving the minimum deposited energy to 60 TeV alleviates much of this tension, since the number of expected 
background tracks only accounts for only ~2 out of the 4 observed events. This can be seen in Fig. 8 where we depict 
the flavor contours for the 60 TeV—3 PeV (left panel) and the 60 TeV—10 PeV (right panel) cases obtained with the 
6 P analyses. As can also be seen in Table I, extending the analysis above 3 PeV yields a dramatic shift in the location 
of the best fit: rather than preferring a dominant Ve component, a strong Vt component becomes favored. As we will 
return to in Sec. V C, this can be understood by the increasing effect of the Glashow resonance as we consider energies 
around ~6.3 PeV. All in all, the low statistics in this range, only 20 observed events (16 showers and 4 tracks), does 
not allow us to reach any strong conclusion regarding the flavor composition and the entire parameter space for the 
standard scenario with averaged oscillations during propagation from the sources (the blue sliver) is allowed within 
la C.L. 

Finally, let us note that the flavor composition remains almost identical to that obtained with the 6 P analyses when 
performing a 7P fit, i.e., when the prompt atmospheric neutrino flux is included with a prior using IceCube limits 
and an additional prior on the atmospheric muon background is implemented. 


C. The Glashow resonance and broken power laws 

In Fig. 8, by considering the cases with i?min = 60 TeV, we show the effect of the upper deposited energy limit 
of our analysis region. By extending this limit from 3 TeV (left panel) to 10 PeV (right panel), the best fit shifts 
dramatically, from (og : ct/i ■ ar)e ~ (1 : 0 : 0)0 to ~ (0 : 0 : 1 ) 0 . This shift also occurs in the > 28 TeV case — albeit 
less dramatically, as the lower-energy events remove some statistical weight from high energies. This difference in the 
best fit (and in the spectral index) is mainly due to the absence of events observed near the Glashow resonance peak. 
For the IceCube best fit assuming the canonical (1:1: 1)0 one expects about 1.2 events above 3 PeV. On the other 
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FIG. 9. Effect of adding a break of one unit in the spectral index of the astrophysical neutrino spectrum at Ei, = 
1 Pe V. Ternary plots of the profile likelihood exclusions of the neutrino flavor composition of the 20 high-energy events seen 
at IceCube in the energy range 60 TeV—10 PeV with no break in the spectrum (left panel) and with a break of Jy = 1 at 
Sj/ = 1 PeV (right panel). We perform 4P analyses: atmospheric backgrounds are fixed as in the 3P case, but we allow the 
spectral index 7 to vary in each case. Same format as Fig. 7. Table I shows the best-fit values for these cases, as well as for 
the 60 TeV—3 PeV energy range. 


hand, for the best fit in the 6 P case for the energy range 60 TeV—10 PeV (almost a purely i^r flux), only ^0.3 events 
are expected above 3 PeV (see bottom right panel of Fig. 5). We further tested this by performing MultiNest runs 
with interactions with electrons artificially removed. This leads to a return of the best-fit region toward (1:0: 0 ) 9 . 

The importance of the lack of events around 6.3 PeV has already been discussed in the literature [122]. It could be 
an indication of a broken power law [59] or even a cutoff [37, 186] in the astrophysical neutrino spectrum. We show 
the effect of adding a break in the power-law spectrum, by setting it at 1 PeV and by modifying the spectral index 
by one unit for higher energies, i.e., such that dy = 'yiEi, > 1 PeV) — "f{Ei, < IPeV) = 1. We perform 4P analyses, 
wherein we fix the atmospheric fluxes, but allow 7 to vary. 

Adding a break of one unit in the spectral index dramatically reduces the expected event rate around the Glashow 
resonance, bringing the best fit closer in line with expectation. It furthermore allows for a harder spectrum in the 
region where the events were observed: when a break is added, the best-fit j{E^, < 1 PeV) goes from 2.48 to 2.34 
in the 60 TeV—3 PeV energy range (2.49 to 2.43 for 60 TeV—10 PeV). The most striking effect is on the flavor 
composition, though. In Fig. 9 we show this for the two cases with a minimum deposited energy of 60 TeV and by 
performing 4P fits. When the energy range is extended up to 10 PeV, adding a break moves the best-fit location 
from (0 : 0.02 : 0.98)0 to (0.75 : 0.25 : 0 ) 0 , as a higher De flux becomes allowed due to the flux suppression above 
E,y = 1 PeV. This shows that if such a break does indeed exist, omitting it when fitting the high-energy events can 
lead to a major mischaracterization of the flavor composition. However, with the current statistics, performing a full 
spectral analysis including the two extra parameters to describe the broken power law (the position of the break and 
the spectral index at energies above it) would be somewhat fruitless. 


D. Track/shower misidentiflcation 


We also consider the possibility that some fraction of the tracks produced in the detector could be misidentified as 
showers. This could occur for instance when a muon neutrino produces a muon close to the borders of the detector, 
but the muon escapes undetected and only the hadronic shower is seen. Although this has been mentioned in a recent 
IceCube analysis with a low-energy threshold at 1 TeV [185], we could reproduce the numbers of the two- and three- 
year analyses above ~30 TeV without including this information [2, 3]. Let us note that this was one of the potential 
explanations to alleviate the tension for the ( 1 : 1 : 1)0 case that was put forward in our previous works [ 120 , 121 ]. 

Here, we show the effect of the misidentiflcation of muon tracks as showers, for both incoming astrophysical and 
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FIG. 10. Effect of including muon track misidentification as shower events. Ternary plots of the profile likelihood 
exclusions of the neutrino flavor composition of the 20 high-energy events seen at IceCube in the energy range 60 TeV—10 PeV 
with a 20% fraction of muon track misidentihcation (left panel) and with a 30% misclassification fraction (right panel). We 
perform 4P analyses. This shows that a less extreme flavor ratio is easily obtained if some of the muon tracks are missed and 
the associated events are identihed as showers. These should be compared with the left panel of Fig. 9, which shows the same 
case but with no misidentihcation of tracks. Same format as Fig. 7. In Table I, we also indicate the best hts for the other 
energy intervals for the case of 20% track misidentihcation. 


atmospheric neutrino events®. We illustrate the effect in the energy range 60 TeV—10 PeV in Fig. 10 by performing 4P 
fits (leaving the spectral index free), which should be compared with the left panel of Fig. 9, where no misidentification 
of tracks was included. We depict the ternary plots of the profile likelihoods in flavor space for the case with 20% 
track/shower misclassification (left panel), and in the case where this fraction reaches 30%® (right panel). The shift 
toward a larger muon neutrino component in the astrophysical flux is clear. In this case, a pure Ve flux, which is 
very close to the best fit when the considered energy range is 60 TeV—3 PeV and if no track misidentification takes 
place, is disfavored at more than Icr C.L. As can also be seen in Table I, where we provide the best fit for other 
analyses, track/shower misidentification allows for a larger astrophysical muon neutrino component, as the expected 
tracks from atmospheric events no longer necessarily swamp the signal events which, on the other hand, consist of 
a larger relative number of showers. And for the 4P analysis, for 20% (30%) track misclassification as showers, this 
results in a 0.82 (0.84) p value for (1:1: 1)0 when the energy interval is 60 TeV—10 PeV (as compared to 0.69 when 
no misidentification occurs). Let us also note that track misidentification does not change the best fit for the spectral 
index, however. 


E. Atmospheric backgrounds 

In our 6P fits, we show what happens when the atmospheric background event rates are allowed to freely float. In 
Fig. 11 we show the one-dimensional profile likelihoods of the background atmospheric neutrinos and muons when all 
six parameters are allowed to vary freely for the four energy intervals under discussion. The best-fit values, posterior 
means and errors are given in Table I. In most cases, a better fit can be obtained with a lower atmospheric component 
than that expected by the IceCube Collaboration, although the Icr errors include the expected values in every case. 
The fit for the 60 TeV—3 PeV interval results in a particularly large number of atmospheric neutrino events. On the 
other hand, a large number of atmospheric muons is predicted in the 60 TeV—10 PeV interval. In both cases, the 
required fluxes would swamp the signal below 60 TeV. This is an indication of the low statistical power of the data 
sample to determine the number of background events and of the fact that we do not use the angular information. 


® The reverse, i.e., a shower being misclassified as a track, is very rare [185]. We also note that a 10% misidentification fraction for 
atmospheric muons is already included in all our fits. 

® After the first version of this work, the IceCube Collaboration released a flavor study that quoted the fraction of track misidentification 
as 30% [187]. Note that differences with respect to the case of 20% track misidentification, which was presented in our work before 
the IceCube preprint appeared on the arXiv, are very small, and IceCube results agree with our findings. Another recent analysis also 
appeared after our first version [188], and it also agrees with our main conclusions. 
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N, N, 


FIG. 11. The p values for the total number of events induced by conventional atmospheric neutrinos (left panel) 
and that of veto-passing atmospheric muons (right panel) in each of the energy ranges considered in this work. They 
are obtained after performing 6P fits. 





FIG. 12. The p values for the number of conventional atmospheric neutrinos (left panel), atmospheric muons 
(middle) and prompt atmospheric neutrinos from charmed meson decays (right panel), for the 7P analyses, i.e., 
when a prompt component is added to our fit and priors on and Np are added as discussed in Sec. IV. 


We then consider the case when a prompt atmospheric neutrino component is added to the analysis, corresponding 
to the 7P case, wherein and Np are subject to priors based on independent measurements. Without this prior, 
the data show a preference for a large prompt component over the astrophysical signal. However, in this work we are 
not including the directional information of the events. Using the veto probability discussed in Sec. Ill, the rate of 
prompt neutrinos from the Southern Hemisphere must be approximately half of that from the Northern Hemisphere, 
and no evidence of such a suppression is supported by data [3]. Therefore, the prior on Np is not only well justified 
but required. 

In Fig. 12, profile likelihoods for the total number of events of conventional atmospheric neutrino (left panel), 
atmospheric muon (middle panel) and prompt atmospheric neutrino (right panel) backgrounds in the 7P case are 
shown. The best-fit points, Bayesian posterior means, and uncertainties are given in Table H. In most cases, a very 
small component from charmed meson decays is favored: fewer than one event in the three-year data set, although 
the IfT C.L. errors can accommodate up to four extra prompt events. Whereas the number of atmospheric muons 
slightly changes as a consequence of the prior imposed on it, the number of events induced by conventional neutrinos 
remains statistically the same as that obtained with the 6P fits. The main effect of an added prompt component on 
other model parameters is to allow for a slightly harder astrophysical component, with a slightly smaller number of 
astrophysical events. All in all, 6P and 7P analyses provide very similar results for the rest of the parameters. 
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VI. CONCLUSIONS 

The first evidence for a high-energy neutrino flux of extraterrestrial origin has been obtained after three years of 
observation by the IceCube experiment [1-3] . In this period of time, a total of 36 events (plus one event whose energy 
and direction cannot be reconstructed) have been detected with deposited energies between ~30 TeV and ~2 PeV. In 
Refs. [120, 121], we showed that by looking only at the flavor composition within a single energy bin (28 TeV—3 PeV), 
there is an apparent tension between the canonical astrophysical signal (1:1: 1)0 and the observed ratio of tracks 
to showers, leading to a best fit for the signal at (1 : 0 : 0 ) 0 . In spite of the relatively small number of observed 
events, interesting conclusions can be drawn from a spectral analysis of the observations. Spectral analyses of the full 
data have been performed, either fixing the flavor composition to (1 : 1 : 1)0 [2, 3, 124, 125] or varying the flavor 
composition but with fixed backgrounds [123]. In this work we performed, for the first time, a detailed spectral analysis 
of the 36 high-energy neutrino events detected by IceCube, where the astrophysical flavor composition, spectral index, 
and normalization, along with the number of background events were left free to vary. Our results are summarized in 
Tables I and II, along with a series of figures, which drive the main discussion of the spectral shape of the spectrum 
as well as of the flavor composition of the astrophysical flux. 

We have shown that the picture includes many more subtleties than what a single-energy bin analysis could reveal. 
By adding spectral information, we have shown that an astrophysical or component may be required for a good fit 
to the higher-energy tracks that were observed. Indeed, these tracks are not well explained by the expected atmospheric 
muon and neutrino fluxes due to their steep spectra. The current data prefer a slightly larger astrophysical — and 
thus smaller atmospheric — component, although there are still too few events for this to be statistically significant. 

Concurrently, we have confirmed that the lack of observed events above 2 PeV carries important information about 
the Ve astrophysical spectrum in the form of missing events around the Glashow resonance (see also Ref. [122]). 
We have done so by means of comparing the results for the energy intervals with a maximum deposited energy of 
3 PeV and those with maximum deposited energy of 10 PeV. Moreover, we have also considered the possibility of 
a break in the astrophysical spectrum and how this would impact our results. In general, considering the lack of 
events above 2 PeV for unbroken power-law spectra implies both a larger Vt flux and a steeper spectrum. This is 
more pronounced in the case of a higher minimum energy, which leaves background events with higher statistical 
power out of the analysis. Nevertheless, if these trends become significant with future data, we are left with a puzzle: 
a preference for a large electron neutrino component between 28 TeV and a few PeV due to a paucity of tracks, 
combined with a lack of electron antineutrinos at high energies, reflected by a dehcit of events around the Glashow 
resonance. The latter could have a simple explanation, such as a break in the power law around 1 PeV, which could 
have a connection with the ultrahigh-energy cosmic-ray paradigm [13, 21, 37, 41]. This would not be surprising, since 
the dominant production mechanism of high-energy neutrinos is thought to be spallation of high-energy cosmic rays 
by photons or protons in the intergalactic medium, yielding neutrinos that are ~10 times less energetic. A different 
dominant production mechanism could also be partly responsible. For example, py-only scattering vastly reduces 
the expected Pe flux, and tt^ decays in dense environments could furthermore lead to a suppressed electron neutrino 
component, as the resulting muon can have time to lose a significant fraction of its energy to the source environment 
before decaying [122]. Conversely, a principal origin in the form of neutron decays leads to a higher Vf. component, 
compounding the problem. Such scenarios have been invoked to explain the lack of tracks in the observed data and 
could plausibly dominate [47]. This tension could be reduced by a two-component astrophysical flux [46], though. At 
present all these possibilities are still allowed by observations. 

Although at present it is not statistically significant, if the best-fit point in flavor space continues to lie outside 
of the source triangle (blue sliver in our ternary plots) where flavor compositions allowed by averaged neutrino 
oscillations lie, we could be confronted with an even bigger puzzle. In principle, the lack of expected tracks could 
suggest either a misunderstanding of the atmospheric backgrounds or a misidentification of tracks as showers. Indeed, 
in Fig. 10 we have shown that a misclassification of tracks as showers would imply a higher likelihood for the source 
triangle and would point to a standard physics origin. More compellingly, this signal could be a hint of exotic physics 
such as in some scenarios of neutrino decay [100, 101, 112-114], pseudo-Dirac neutrinos [104, 109, 113, 115], sterile 
neutrino altered dispersion relations due to shortcuts in an extra dimension [70], TeV gravity [71], T violation [58], 
scenarios with a color octet neutrino [57, 69], Lorentz violation [51, 59, 63, 75], leptoquarks [53] or neutrino-secret 
interactions [60, 61, 64, 67, 68, 73]. 

The era of neutrino astronomy has just started, triggered by the first ever detected high-energy neutrinos in IceCube. 
After three years of data, a rate of 10 — 20 events per year is not high. However, even with such a small sample, 
compelling questions are raised on the production, propagation and detection of neutrinos at high energies. More 
events are required: even though the next data release might shed some light on some of the questions elicited 
by our analyses, we see an important need for experiments with higher sensitivity and exposure. In this sense a 
future high-energy extension of the IceCube detector [77] and the planned KM3NeT telescope [189] will be of crucial 
importance. 
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Appendix A: Neutrino cross sections 


For the neutrino-nucleon and antineutrino-nucleon differential cross sections we use the NUSIGMA Monte Carlo 
code [190], which uses the CTEQ 6 -DIS parton distribution functions [191, 192]. The NC and CC differential cross 
sections for neutrinos (and similarly for antineutrinos) of a given flavor off a target molecule with mass, atomic 
and neutron numbers A, Z, and N^ respectively, is 
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(Al) 


where A = 18, Z = 10 and N = 8 correspond to water (for an isoscalar target, A/2 = Z = N) and /dy are 

the NC/CC neutrino-proton(neutron) cross sections. 

For the case of neutrino interactions with electrons per molecule, one needs to include a factor Z/A (as for the 
neutrino-proton cross section). The relevant cross sections are given by [127, 193] 
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where £ = y or t, me and are the electron and muon masses, s = 2meE^^ Gp is the Fermi constant, gp = 
sin^ 9w — 1/2 and gp = sin^ 9^ are the chiral couplings of the Z boson to the electron, 9w is the weak mixing angle, 
Miv and Mz are the gauge boson masses, Tw = 2.085 GeV is the W boson width, T{W —>■ hadrons) = 0.6760 and 
r(lT —>■ Dfj,y) = 0.1057 [156]. The inelasticity y is given by Eg = y E^, for all flavors, £ = {e,y,T}. Note that for the 
case of 9e interactions via the Glashow resonance with a pure hadronic hnal state, the true deposited energy is equal 
to the neutrino energy, so only the total cross section is relevant, see Eq. (B28). 
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Appendix B: Full differential spectra formulas 

For an astrophysical flux /dE^, (assuming the same spectrum and normalization for neutrinos and antineutrinos) 
with flavor combination at Earth {ag : : Ot}®, the differential event spectra for topology k = {sh, tr} is given 

by" 


dNk,a 

^-^dep 




dNg’‘^ 

^-^dep 


(Bl) 


In this appendix we present the full expressions for the differential spectra, dN^’^/dE^gp, for showers {k = sh) and 
tracks {k = tr) for each channel £ = {e, /i, r}, and for an incoming neutrino flux of type / (astrophysical, conventional 
atmospheric or prompt atmospheric flux, / = {a, z/, p}). The final spectra result from the sum of all the partial 
contributions from the different processes and are given by 
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dN^h,NCJ 

dEfigp 


dE, 


dN; 


dep 

sh,NC,/ 


dE, 


dN; 


dep 

sh,NC,/ 


dE, 


dep 


dN, 


■tr,e./ 


+ 


dN, 


tr,e./ 


^T/dep 

dNl^,cc,f 

^-^dep 

dA,)r.CC./ 


dE, 


dep 


dN, 


tr.CCJ 


^-^dep 

^^tr.CC,/ 


dE, 


dep 


dE, 


dep 


dN^hfiCJ rfTV^h.CC./ ^^sh.e./ 

dEg^gp dEg^gp dEg^gp 

_ |_ 

dE(^Q^ dE(\Qp 

^jYsh,cc,/ 

u-r ^ u-r u-r 

^^dep dE^QT^ ^^dep 


^^dep 

^^dep 


dNp^d 

__ 

^-^dep 


dNp^d 

dEg\gp 


dN, 


sh,e, 


dE, 


dep 


/ 

- + 


*^-^dep 




dE, 


dep 


,(B2) 

(B3) 

(B4) 

(B5) 

(B6) 

(B7) 


and all the partial contributions from the different processes are detailed below in this appendix. 

Showers are induced by both Vg and CC interactions with nucleons, as well as by NC interactions of neutrinos 
of all three flavors with nucleons. The differential shower spectrum, in terms of the measured deposited energy, by 
NC interactions for neutrinos (and analogously for antineutrinos) with flavor £ reads 


^^sh.NCJ 

dEgigp 


= TNa 



dEgg Attl^{E,,) 


d4>UE,) 

dEy 


[ dyMge{E^^)R{E^^,Eg,gp,aiE^^))E£ldEyl 

Jo 


(B8) 


where T = 988 days is the time of data taking, Na = 6.022 x 10“^^ g“" and E^^y = {Eg, — E^) is the hadronic shower 
energy, with E'g, the energy of the outgoing neutrino. The attenuation and regeneration factor due to the absorption of 
neutrinos when traversing the Earth is given by Attl^{E,g) (see Sec. II A) and the detector effective mass as a function 
of the true deposited energy by Meff(Etrue) (see Sec. IIC). The incoming neutrino flux of type / (/ = {a, v,p} for the 
astrophysical, atmospheric neutrino and atmospheric muon flux, respectively) is d(f)l^{Eg,)/dE,g, and da^^/dy is the 
neutrino-nucleon NC differential cross section. The energy resolution function is given by R{EtTue, E,igp, a{Etgue))- 
We assume the uncertainty on the true EM-equivalent deposited energy, Etme, to be given by the error on the 
measured EM-equivalent deposited energy, Edep, and we perform two single-parameter (e) fits, for Edep < Etme and 
Adep > A'true; with the function a{Etgg,g) = eEtrue within the observed energy range using the 36 (shower and track) 
events detected by IceCube after 988 days. Since upper and lower errors are different, we represent R with two 
half-Gaussians, 


R{EtTue: Egigp, (j{Ef,ggig'^^ — 


(Etrue-Edep)^ 


V2 


TTCT^ 


(B9) 
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Note that for the conventional and prompt atmospheric neutrino fluxes the relative contributions from each flavor, and for neutrinos 
and antineutrinos, are obtained by using the corresponding flux for each case. 
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where the best fits for the lower and upper dispersions are 




0.121 E^Yue If -^dep ^ ^true 

0.125 £^true if -^dep ^ -^true ■ 


(BIO) 


The differential shower spectrum produced after CC neutrino-nucleon interactions has contributions from Ve and 
Vr (and similarly from Dg and Pt)- In the case of Ve (and analogously Pe) CC interactions, the differential shower 
spectrum in terms of the measured deposited energy is given by 


dilldep 


TNa 



dEi, Attl^{E„) 


d(l>lSEv) 

dE^ 




CC 


rCC\ 


dc^u^[E^,y) 

dy 


(Bll) 


In the case of Vr (and analogously i>t) CC interactions, the differential shower spectrum in terms of the measured 
deposited energy reads 


diVfh’CC./ 


dE, 


dep 


= TNa dE.Attl (E,) 


dfpl^ {E^) 


dEi, 


lo dy 


dz 


(B12) 


X ^ {DAEr)M,s{E^;^)R{E^;^,Ei,p,a{E^;^)) + il-DAEr))M,s{Eh)R{Eh,E^,p,a{Eh))) 


c=h,e 


where dnc{z)/dz is the energy distribution of the daughter Vr or e with energy [z = E^^^t/Et) or E^^t (z = 

Ee^r/Er) from T decay via the hadronic or electronic channel (c = {h, e}), respectively [140]. The fraction of tau 
leptons decaying inside the detector, Dt(Et), is dehned in Eq. (27) and a ht in terms of a Fade approximant, with 
an accuracy at the 0.1% level, is given by 


Dr{E,) = 


1+Pi {Er/lOPeY) 


1 + qi {Er/10PeV) + q 2 {Er/lOPeY) 


2 ’ 


(B13) 


with Pi = 0.883, qi = 1.66 and q 2 = 1.15. 

On the other hand, tracks are produced in CC 12 ^ and I'r interactions (followed by the r decay into For 

12 ^ (and analogously P^) CC interactions, the differential track spectrum in terms of the measured deposited energy 
reads 


_ 

^^dep 


TNa 



dE^ Attl^{E^) 


d(t>l^{E^) 

dEi, 



dy Meff(E^'^) RiE^^, Edep, AE^^)) 


dcr^^{E,,,y) 

dy 


(B14) 


The differential track spectrum from v-r (and analogously from v^) CC interactions is 


^^tr.CC,/ 


dE, 


dep 


= TNa 


dE^Attl {El,) 


d<Pt^{E^) 


dEi, 


C dy 

lo dy 


(B15) 


X ^ dz {D,{Er)M,s{E^^)R{E^^,E,i,p,a{E^^)) + {l-DriE,))M,g{Eh)R{Eh,Ed,p,a{Eh))) 

where dnfj_{z)/dz is the energy distribution of the daughter p with energy {z = E^j^^r/Er) from r decay [140]. 
The true EM-equivalent deposited energies are given by 


Eh = Eh{E„y) E„y , 

Ee = E^{l-y) , (with £ = {e, p, t}) , 

E^^ = Eh , 

Ee^ = Eh + Eg , 

E^^ = Eh + Ffi {Efj_ + a/b) , 

=Eh + Fh{E^ (1 - z)) E^ (1 - z) , 

Et!^ = Eh + Et z , 

E?lll = Eh + F^^t{Et) {Et z + a/b) , 

where Fh{Ex) is defined in Eq. (20), a and b in Eq. (21), F^ = 0.119 in Eq. (25) and F^^t{Et) in Eq. (28). 


(B16) 

(B17) 

(B18) 

(B19) 

(B20) 

(B21) 

(B22) 

(B23) 
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The function Fh{Ex) can be parametrized as [194] 


= , (B24) 

where the values of the parameters resulting from a ht to simulations of hadronic showers induced by NC neutrino- 
nucleon interactions are /o = 0.467, E^ = 0.399 GeV and m = 0.130 [155]. 

As for a and 6, we use tabulated data for the muon loss rate in ice [156] and perform a fit taking both parameters 
as constants. The resulting values are: a = 0.206 GeV/m and b = 3.21 x 10“"^ m“^. This ht is accurate at the few 
percent level above 400 GeV and below the percent level above 40 TeV. 

The average fraction of energy lost along a track of a muon produced in a tau decay after a v-r or GC interaction 
inside the detector, E^^t^Et), is dehned in Eq. (28) and a ht in terms of a Fade approximant, with an accuracy at 
the 0.2% level, is given by 


Ffj,.r{Er) = F^ 


_ 1+pi (Er/lOPeY) _ 

1 + qi (E^/lOPeV) -f 92 {Er/lOPeVf ’ 


(B25) 


with Pi = 0.984, qi = 1.01 and 92 = 1-03. 

The contributions to the differential shower spectrum from neutrino or antineutrino interactions with electrons from 
diherent processes are given by 




dEdep „ „ 

where Vi = {ve, Vfi, Vt, Ve, S't}, and 


v Jo 


dy 


dNl 


sh,e,f 


dE, 


dep 


■TNa [ dE,AuiiE,)^^y^ [ 

Jo Jo dy 


c.— h,e 


(B26) 


(B27) 


dN': 


sh,ej 


Ve,h 


dE, 


dep 


TNa dE, AttiSE.) M,s{El^^) R{El^, Edep, al^„{E,) , (B28) 


^-^dep 


'^=TNa r dE,Attl{E,)^^^^ 
ep Jo dEi, Jq dy 

JO 


dnc{z) 

dz 


(B29) 


The contributions to the differential track spectrum from neutrino or antineutrino interactions with electrons from 
diherent processes are given by 


dNl\Y 


dE, 


dep 


r ^ (]t) (z\ 


where ui = {vs^Vt}, and 


(B30) 


=TNa dE, Atti^iE,) dyM,s{E;) R{E;, Edep.., a{Ef,)) 


(B31) 


with VI = {Vf,, Vft}. 
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Deposited EM-Equivalent Energy in Detector [TeV] 


FIG. 13. Event spectra in IceCube, as a function of the EM-equivalent deposited energy, of tracks and showers 
for each flavor after 988 days (summing neutrinos and antineutrinos), for an isotropic power-law spectrum, per flavor, 
d(f>°'/dE„ = 1.5 X 10“® (Ujy/lOO TeV)“°'® GeVcm“^ s“^ sr~^ [3]. Showers (tracks) are depicted as solid (dashed) histograms. 
The contributions from i>e + he, Ufj, + and Vr + Vr are represented by red, green and blue histograms, respectively. 


The true EM-equivalent deposited energies for all the neutrino-electron processes are 


El = E,y , 

(B32) 

E-,^^ = Fh{E,y{l-z))E,y{l-z) , 

(B33) 

Er,e = E,,yz , 

(B34) 

El^ = Fn{E,)E,, , 

(B35) 

Er.f, = F^^riE,,y){E,,yz + a/b) , 

(B36) 

El = F^{E,y + a/h) . 

(B37) 


In Fig. 13 we show the event spectra of showers and tracks for each flavor (summing neutrinos and antineutrinos) for 
the best ht IceCube spectra [3], ie., jdEy = 1.5 x 10“® (E^/100TeV)“°-^ GeVcm“^ s“^ sr“^, per flavor. The 

effect of the Glashow resonance on the f'e-induced event spectra is clearly visible in the red histograms. The shower 
spectrum for Vf, -|- (red solid histogram) shows a bump above a few PeV and the resonant interactions of i>(, []but 
also the no-resonant interactions of Ve\ see Eq. (A6)] with electrons also give rise to tracks (red dashed histogram), 
via the the leptonic decay of the produced W bosons. We also note the similar shape of all the event distributions 
(except from the red dashed histogram), as a function of the EM-equivalent deposited energy, below a few PeV. 
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